Efficient evaluation of partition functions of frustrated and inhomogeneous spin 
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We present a numerical method to evaluate partition functions and associated correlation func- 
tions of inhomogeneous 2-D classical spin systems and 1-D quantum spin systems. The method is 
scalable and has a controlled error. We illustrate the algorithm by calculating the finite-temperature 
properties of bosonic particles in 1-D optical lattices, as realized in current experiments. 
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The study of thermodynamic properties of 2-D classi- 
cal spin systems has a very long and fruitful history. On 
the one hand, it has provided valuable insights in the the- 
ory of magnetism and phase transitions. On the other, it 
has allowed us to describe 1-D quantum systems at finite 
temperature by virtue of the Suzuki-Trotter decomposi- 
tion. 

As very few spin models are exactly solvable, many dif- 
ferent approximate methods have been proposed to cal- 
culate the associated partition functions. Monte Carlo 
seems to be the method of choice for non-frustrated 
systems, but fails in the description of frustrated and 
fermionic systems due to the notorious sign problem 0. 
No such sign problem occurs in the method developed by 
Nishino |2] , where the largest eigenvalue of the transfer 
matrix of the classical spin model can be approximated 
by using a variation of the density matrix renormalization 
group (DMRG) approach 0, By making use of the 
Suzuki-Trotter decomposition, this method has also been 
used to calculate the free energy of translational invari- 
ant 1-D quantum systems 0, |a ■ The main restriction 
of this method is that it cannot be applied in situations 
in which the number of particles is finite and/or the sys- 
tem is not homogeneous. The method may also become 
ill conditioned when the transfer matrix is not hermitian. 
Finally, in the case of 1-D quantum systems, a recent de- 
velopment 1^ allows us to overcome these problems by 
extending the concept of matrix product states [ifl, UM 
to mixed states, e.g. by using the idea of purification of 
states 0. This method is, however, specially designed 
for 1-D quantum systems, and cannot be extended to 
classical 2-D models. 

Here we take a completely different approach which 
allows us to overcome the drawbacks of the above men- 
tioned methods in both 2-D classical and 1-D quantum 
systems. We achieve this by evaluating the associated 
partition and correlation functions directly. The main 
advantages of this method are that the approximations 
made are very well controlled, that it applies to frus- 
trated, inhomogeneous and finite classical and quantum 
systems, and that it can be generalized to higher dimen- 
sions. We will illustrate its performance with the system 



of strongly interacting bosonic atoms in optical lattices, 
a problem which has attracted a lot of interest in the 
atomic physics community in the last few years due to 
the recent experimental achievements [l^ [iSl [ij] . In this 
system, 20-100 atoms are trapped by the combination of 
a periodic and a harmonic (i.e. inhomogeneous) poten- 
tial created by lasers in 1-D and at a finite temperature. 
In the so-called Tonks-Girardeau limit the prob- 
lem can be exactly solved via fermionization [l^ . and 
thus this provides us with a reliable benchmark for our 
method. Outside this limit, we are able to reproduce 
certain features experimentally observed |l2j . 

Our method relies in reexpressing the partition and 
correlation functions as a contraction of a collection of 
4-index tensors, which are disposed according to a 2-D 
configuration. We will perform this task for both 2-D 
classical and 1-D quantum systems. We will then show 
how the methods introduced in can be used to ap- 
proximate these contractions in a controlled way, and 
thus lead to a scalable algorithm for the evaluation of 
the quantities of interest. 

1.- 2-D classical systems: Let us consider first the par- 
tition function of an inhomogeneous classical 2-D n-level 
spin system on a Li x L2 lattice. For simplicity we will 
concentrate on a square and nearest-neighbor interac- 
tions, although our method can be easily extended to 
other short-range situations. We have 

Z= J2 exp[-/3i7(x",...,a;^^^^)] , 

a;ii,...,a;'^i^2 



where 



is the Hamiltonian, x^^ = 1, . . . , n and [3 is the inverse 
temperature. The singular value decomposition allows 
us to write 

n 

exp [-l3Hl^{x,y)\ = ^ /^^Ja;)5^i(y), 

a=l 
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with q € {i, — >}. Defining the tensors 

n 



^Irud 



x = l 



the partition function can now be calculated by contract- 
ing all 4-index tensors X^^ arranged on a square lattice 
in such a way that, e.g., the indices l,r,u,d of X*-' are 
contracted with the indices r,l,d,u oi the respective ten- 
sors X^'-'-i, X^'-'+i, X^+i'-J'. In order to determine 
the expectation value of a general operator of the form 
0{{x^^}) = ZjJ-j 0*-'(x*-'), one just has to replace each 
tensor X*-' by 



x=l 



2.- 1-D quantum systems: We consider the partition 
function of an inhomogeneous 1-D quantum system com- 
posed of L ri-level systems, 

Z = Trexp(-/3i/). 

It is always possible to write the Hamiltonian H as & 
sum H = Hk with each part consisting of a sum of 
commuting terms. Let us, for simplicity, assume that 
H = Hi + H2 and that only local and 2-body nearest 
neighbor interactions occur, i.e. Hk = J2i ^k^^^ ^^-'^ 
0^''+\ O^'^'+^j = 0, with i,j = 1,...,L. The more 

general case can be treated in a similar way. Let us now 
consider a decomposition 



exp 



^ka 



'T, 



i+l 



(1) 



The singular value decomposition guarantees the exis- 
tence of such an expression with k < n^. As we will see 
later, a smart choice oi H ^ Hk can typically de- 
crease K drastically. Making use of the Suzuki-Trotter 
formula El 



^ = Tr(nexp(-^i7. 
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it can be readily seen that the partition function can 
again be calculated by contracting a collection of 4-index 
tensors X'^^ defined as 



')ud — ^ ll^lr-^2l'^2r' 



[ud] 



where the indices (Z, I') and (r, r') are combined to yield a 
single index that may assume values ranging from 1 to . 
Note that now the tensors X'-* and X^ ^ coincide, and 
that the indices u of the first and d of the last row have 
to be contracted with each other as well, which corre- 
sponds to a classical spin system with periodic boundary 



conditions in the vertical direction. A general expecta- 
tion value of an operator of the form O = ZO^ (g) • ■ • (g) 
can also be reexpressed as a contraction of tensors with 
the same structure: it is merely required to replace each 
tensor X^^ in the first row by 



■^(U'){rr')ud 



ud] 



3.- Tensor contraction: In the following, we adapt the 
algorithm introduced in 17] in order to contract the ten- 
sors X'^^ introduced above in a controlled way. The main 
idea is to express the objects resulting from the contrac- 
tion of tensors along the first and last column in the 2 D 
configuration as matrix product states (MPS) and those 
obtained along the columns 2,3,..., L—1 as matrix prod- 
uct operators (MPO) :,S|. More precisely, we define 

m 
rn 

^ 1 ■ ■ ■ ^ J\-/ — 1 

h,ru... = l 

where m ~ n for 2-D classical systems and m = 
for 1-D quantum systems. These MPS and MPOs are 
associated to a chain of M m-dimensional systems and 
their virtual dimension amounts to D = n. Note that 
for 2-D classical systems the first and last matrices un- 
der the trace in the MPS and MPO reduce to vectors. 
The partition function (and similarly other correlation 
functions) reads Z = ( X^ | X^ • • • X^"^ | X^ ). Evaluat- 
ing this expression iteratively by calculating step by step 
( X^ I ( X^-i I X^ for j 2, . . . , L - 1 fails because the 
virtual dimension of the MPS ( X-* | increases exponen- 
tially with j. A way to circumvent this problem is to 
replace in each iterative step the MPS ( X-* | by a MPS 

( X"* I with a reduced virtual dimension D that approx- 
imates the state (X-* | best in the sense that the norm 

6K := II (X-' I - (X'' III is minimized. Due to the fact 
that this cost function is multiquadratic in the variables 
of the MPS, this minimization can be carried out very 
efficiently p, 0, ^3| ; the exponential increase of the vir- 
tual dimension can hence be prevented and the iterative 
evaluation of Z becomes tractable, such that an approx- 
imation to the partition function can be obtained from 
Z ~ ( X I X^ ) . The accuracy of this approximation 
depends only on the choice of the reduced dimension D 
and the approximation becomes exact for D > D^. As 
the norm SK can be calculated at each step, D can be in- 
creased dynamically if the obtained accuracy is not large 
enough. In the worst case scenario, such as in the NP- 
complete Ising spin glasses 0|, D will probably have to 



3 



grow exponentially in L for a fixed precision of the par- 
tition function. But in less pathological cases it seems 
that D only has to grow polynomially in L; indeed, the 
success of the methods developed by Nishino Q in the 
translational invariant case indicate that even a constant 
D will produce very reliable results. 

4-.- Illustration: Bosons in optical lattices: A system 
of trapped bosonic particles in a 1-D optical lattice of L 
sites is described by the Bose-Hubbard Hamiltonian 



H 



L-l 

-J^{a\ai+i + h.c.' 
i=i 



U 



2 ^ 

i=l 



where a| and ai are the creation and annihilation oper- 
ators on site i and = a|ai is the number operator. 
This Hamiltonian describes the interplay between the ki- 
netic energy due to the next-neighbor hopping with am- 
plitude J and the repulsive on-site interaction U of the 
particles. The last term in the Hamiltonian models the 
harmonic confinement of mag nitudeFj = Voii-io)"^- The 
variation of the ratio U/J drives a phase-transition be- 
tween the Mott-insulating and the superfluid phase, char- 
acterized by localized and delocalized particles respec- 
tively ^0] . Experimentally, the variation oiU / J can be 
realized by tuning the depth of the optical lattice [lflll2ll |. 
On the other hand, one typically measures directly the 
momentum distribution by letting the atomic gas expand 
and then measuring the density distribution. Thus, we 
will be mainly interested here in the (quasi)-momentum 
distribution 



nk 



= 7 E {alas)e^'^' 



{r-s)/L 



r,s — 1 



Our goal is now to study with our numerical method 
the finite-temperature properties of this system for dif- 
ferent ratios U / J . We thereby assume that the system 
is in a thermal state corresponding to a grand canonical 
ensemble with chemical potential /x, such that the par- 
tition function is obtained as, Z — Tr e^'^*^^~^^^ Here, 
N = X^i^i represents the total number of particles. 
For the numerical study, we assume a maximal particle- 
number q per lattice site, such that we can project the 
Hamiltonian H on the subspace spanned by Fock-states 
with particle-numbers per site ranging from to q. The 
projected Hamiltonian H then describes a chain of L 
spins, with each spin acting on a Hilbert-space of dimen- 
sion n ~ q+\. A Trotter decomposition that turned out 
to be advantageous for this case is 



M 



1 

.L-l 



(2) 



with H ^ Hr + Hs + Ht, Hr^ J2i=i R^'^R^'+^\ 
a„ = - a,), r(^) = in,(n, - 1) + V^h^ 
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FIG. 1: Density and (quasi)-momentuni distribution in the 
Tonks-Girardeau gas limit, plotted for j3J = 1, L = 40, 
= 21 and Vo/J = 0.034. The dots (crosses) represent 
the numerical results for D — 2 {D — 8) and the solid line 
illustrates the exact results. From the insets, the error of the 
numerical results can be gathered. 



a|, Oi and Ui 



and V = e~2M^Re~2M^sg-2M(^T-fiN) 

thereby denote the projections of the creation, the an- 
nihilation and the number operators a|, ai and Ui on 
the g-particle subspace. The decomposition ^ of all 
two-particle operators in expression ((SJ then straightfor- 
wardly leads to a set of 4- index tensors Xl^^^^, with in- 
dices I and r ranging from 1 to (g H- 1)^ and indices u 
and d ranging from 1 to g -|- 1. Note that the typical sec- 
ond order Trotter decomposition with H = Hcvcn + ^^odd 
would make the indices I and r range from 1 to (q + 1)^. 

Let us start out by considering the limit U/ J ^ oo 
in which double occupation of single lattice sites is pre- 
vented and the particles in the lattice form a Tonks- 
Girardeau gas In this limit, the Bose-Hubbard 
Hamiltonian maps to the Hamiltonian of the exactly solv- 
able (inhomogeneous) XX-model, which allows to bench- 
mark our algorithm. The comparison of our numerical 
results to the exact results can be gathered from fig. ^ 
Here, the density and the (quasi)-momentum distribu- 
tion are considered for the special case PJ = 1,L = 40, 

= 21 and Vq/J = 0.034. The numerical results shown 
have been obtained for Trotter-number M — 10 and two 
different reduced virtual dimensions D ^ 2 and D ~ 8. 
The norm 5K was of order 10^^ for D = 2 and 10^^ 
ioi D — 8 |26|. From the insets, it can be gathered that 
the error of the numerical calculations is already very 
small for D = 2 (of order 10~^) and decreases signifi- 
cantly ioi D — 8. This error can be decreased further by 
increasing the Trotter- number M . 

As the ratio U/J becomes finite, the system becomes 
physically more interesting, but lacks an exact mathe- 
matical solution. In order to judge the reliability of our 
numerical solutions in this case, we check the convergence 
with respect to the free parameters of our algorithm [q, 
D and M). As an illustration, the convergence with re- 
spect to the parameter q is shown in figure [3 In this 
figure, the density and the (quasi)-momentum distribu- 
tion are plotted for g = 2,3 and 4. We thereby assume 
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FIG. 2: Density and (quasi)-momentum distributions for in- 
teraction strengths U/J — 4 and 8. Here, (3 J = 1, L — 40, 
= 21 and M — 10. Numerical results were obtained 
for q = 2 (plus-signs), g = 3 (crosses) and 5 = 4 (solid line). 
For comparison, the distributions for U/J — (dotted hues) 
and (7/J — > cxD (dash-dotted lines) are also included. 



FIG. 3: FWHM of the (quasi)-momentum distribution as a 
function of U/J, calculated for temperatures /3J = 0.5 (plus- 
sign),/? J — 1 (crosses) and (3 J — 2 (dots). The corresponding 
(quasi)-momentum distributions for U/J — 2 and U/J = 8 
are illustrated in the plots at the right-hand side. 



that f3J = 1, i = 40 and = 21 and consider interaction 
strengths U/J — 4 and 8. The harmonic potential Vq is 
chosen in a way to describe Rb-atoms in a harmonic trap 
of frequency Hz (along the lines of ^3)- note that we 
have taken into account that changes of the ratio U/J are 
obtained from changes in both the on-site interaction U 
and the hopping amplitude J due to variations of the 
depth of the optical lattice. The numerical calculations 
have been performed with M — 10 and D = q + 1. From 
the figure it can be gathered that convergence with re- 
spect to q is achieved for q > 3. 

We now use our numerical algorithm to study a physi- 
cal property of interacting bosons in an optical lattice, 
namely the full width at half maximum (FWHM) of 
the (quasi)-momentum distribution. It has been pre- 
dicted that the FWHM shows a kink at zero tempera- 
ture This kink is an indication for a Mott- 
superfiuid transition, since the FWHM is directly related 
to the inverse correlation length. Experiments have also 
revealed this kink they are, however, performed 
at finite temperature, something we can study with our 
algorithm. In figure 13 we plot the numerical results for 
the FWHM as a function of U/ J for three different (in- 
verse) temperatures PJ = 0.5,1 and 2. The physical 
parameters L, TV and Vq are thereby chosen as in the 
previous case. The numerical results have been obtained 
for M = 10, g = 4 and D — q+1. For each temperature, 
three different regions can be distinguished: the super- 
fluid region with constant FWHM, the Mott-region with 
linearly increasing FWHM and an intermediate region in 
which both phases coexist. The value U/J at which the 
Mott-region starts increases with increasing temperature, 
which is consistent with the criteria U ^ fesT, J for the 
appearance of the Mott-phase. This behaviour could be 



easily observed in present experiments. 

In summary, we have presented a numerical method 
for the investigation of thermal states of inhomogeneous 
2-D classical and 1-D quantum systems. We have illus- 
trated the usefulness of this method by applying it to a 
system of trapped bosonic particles in a 1~D optical lat- 
tice - which is of current experimental interest. For this 
system, we have studied the error and the convergence 
of the method. In addition, we have used the method 
to study some physical properties of this system and we 
have obtained results which can be verified experimen- 
tally 
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